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1.0 Introduction 



■ •'.-A l.'Bkary 

SCf^00L 

■’•lOR'rnrt,./ CA a3243.£;-,01 

The physical processes o-f the earth's atmosphere can be 

modelled by a system o-f hydrodynamic equations. This system of 
equations cannot be solved directly unless many simplifying as- 
sumptions are made, severely limiting how realistically the 
actual atmospheric processes are simulated. In order to produce 
accurate weather forecasts the full system of equations must be 
solved in four dimensions. In practice global weather forecasts 
are produced at various meteorol ogical canters around the world 
by treating these equations as an initial value problem and 
integrating forward in time to produce forecasts. The solution 
of this problem requires the use of advanced vector processors as 
the number of computations involved is staggering. The forecast- 
ing problem was formulated quite succinctly by the Norwegian 
pioneer in weather forecasting, V. Bjerknes C13, when he defined 
the necessary and sufficient conditions for a successful system 
in an article written in 1922 to bes 

i) A sufficiently accurate knowledge of the state of the 
atmosphere at the initial time. 

ii) A sufficiently accurate knowledge of the laws according 
to which one state of the atmosphere develops from another. 



Bjerknes' discoveries of the hydrodynamical nature of the 
weather problem led several European nations, especially the 
Scandinavian countries, to begin collecting observations of the 
state of the atmosphere. This data collection led L. F. 
Richardson C2D to try describing initial conditions from a hand 
analysis and projecting the state of the atmosphere to the future 
from the hydrodynamical equations. The task was monumental as 
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Richardson estimated that a warehouse o-f 64,000 people using the 
mechanical calculator of the day could Just forecast the state of 
the atmosphere at the rate that it was actually evolving. Unfor- 
tunately, many factors, discovered later in the 1940 's, kept 
Richardson from making a successful forecast. 

The magnitude of the weather forecasting problem required 
the development of electronic computers for even simple 
solutions. The Electronic Numerical Integrator and Computer 
(ENIAC) developed in the late 1940's allowed Charney, Fjortoft, 
and Von Neumann C3II to succeed in making a reasonable 24 hour 
forecast. Their hydrodynamical model was simplified to filter 
gravity waves while allowing weather patterns to develop in a 
manner similar to that observed in the atmosphere. Their initial 
conditions of pressure heights were derived by hand and the 
result gridded and typed into the computer. 

With the rapid development of computers over the past thirty 
years, it has become possible to use numerical techniques to 
integrate the full set of hydrodynamic equations forward in time 
to produce improved weather forecasts. As V. Bjerknes predicted, 
accurate forecasts require more than Just accurate treatment of 
the physical processes of the atmosphere, they also require 
accurate specifications of the initial state from non-uniformly 
located observations. Panofsky C4J, Bergthorsson and Odds C5J, 
and Cressman C61 pioneered methods to use the computer to obtain 
a weather analysis from observational data. This process of 
combining observation values with a background field is called 
objective analysis. For the most part, these original objective 
analysis techniques were weighted average schemes that depended 
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upon proper specification of several parameters, usually obtained 
in an ad hoc way. Today, most of the world's weather centers use 
statistical abjective analysis techniques based on the work of 
Gandin C73 to provide initial conditions for their atmospheric 
forecast models. 

In practice two sources of information are combined to 
produce an objective analysis: observations of atmospheric 
variables and a forecast made by the atmospheric model from a 
previous analysis. The forecast is commonly called the 'first 
guess' to the analysis or the 'background'. Because the forecast 
is hardly a guess, the term 'forecast background' is used in this 
paper to emphasize that the background to the analysis was 
derived from a forecast made earlier. The observations of 
temperature, wind, and moisture are made by in situ instruments 
attached to balloons, aircraft, and ships or from remote instru- 
ments aboard satellites or on the earth's surface. The result is 
a collection of observations of varying degrees of accuracy taken 
at various times. The statistical analysis schemes have been 
designed to 'optimally' combine these observations with the fore- 
cast background to produce the initial conditions required by the 
numerical forecast model. The optimality of these schemes 
directly depends upon how well the statistical properties of the 
errors of the forecast background and the observations are 
defined. In practice, the schemes are multivariate in the sense 
that they are used to simultaneously analyze multiple related 
dependent variables from measured values. 

In this paper we will deal with the representation of the 
statistics of the forecast background error. In particular, the 
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modelling of the spatial autocovariance of the error for the 
primary variable is examined. Early versions of Gandin's method 
used a simple exponential function to model the autocovariance. 
Although this model is simple, it failed to be sufficiently 
flexible to describe details of the statistics derived from 
actual data. A search of the literature revealed that many 
models are available, but tests of their abilities in fitting 
background statistics for an actual forecast model have not been 
conducted. The mathematical and precision limitations of various 
models have been determined and are described in this paper. 

Optimum interpolation (01), which is sometimes more prac- 
tically referred to as statistical interpolation <SI), is applied 
to compute the corrections to the background field. This is done 
by first interpol ating the background to the non-unif ormiy lo- 
cated observation locations, and then computing the difference 
between the observed and background value. If observations were 
exact, this would be the background error measured at discrete 
locations. These measurements of background error are then used 
by 01 to compute a correction field on a uniform grid, which is 
added to the background to produce the analysis. 

The full development of the equations for a multivariate 
application of 01 is given in several papers, including 
Rutherford C81 , Schlatter C91, Schlatter, et. al . ClOl, Bergman 
Cl 13, Lorenc C123, Thi^baux C133, and Thi^baux and Redder C143. 

A brief outline of the method is given in the following. For a 
collection of estimates of the error at scattered points, it is 
desired to .estimate the value of the error at the grid points. 

□I approximates these values by a linear combination of the known 
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values de-fined so that the expected mean squared error over some 
ensemble oHF realizations is minimized. This requires that the 
statistical properties (covariances between variables) be known. 
Stationarity (independence ot the particular grid point) o-f the 
statistical parameters is required tor a tractable problem. The 
weights used in the linear combination are obtained trom the 
solution ot a certain system ot linear equations, the coetticient 
matrix being the matrix ot covariances between the background 
plus observed errors at the observation points. The positive 
detiniteness ot this matrix plays an important role, both the- 
oretically and computationally. 

A discussion ot multivariate covariance tunctions, proper- 
ties they must satisty, and methods ot obtaining such tunctions 
are discussed in section 2. Experiments with several tamilies ot 
covariance tunctions in titting background error statistics and 
the resulting pertormance in a statistical interpolation scheme 
are described in section 3. Section 4 summarizes the results 
and suggests some tuture work. 

2.0 Multivariate Covariance Functions 

2.1 General Development 

The theory ot covariance tunctions and that ot positive 
detinite tunctions go hand-in-hand. Positive detiniteness ot 
matrices such as occur in our application are equivalent to the 
spatial covariance tunction tor the background errors being posi- 
tive detinite. Positive detinite tunctions are characteri zed by 
Bochner's Theorem C151, which states that a tunction is positive 
(semi ) det ini te it and' only it its Fourier transform is. 
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nonnegative. Alternatively, the covariance -function is the 
Fourier (cosine) transform of a probability density (nonnegative) 
function. Because of the application, our interest is in posi- 
tive definite functions that are smooth in the sense that certain 
partial derivatives exist. An excellent reference for positive 
definite functions is Stewart C161. 

For completeness, a derivation of the covariance functions 
for variables related through differentiation is given here. 
Suppose that it is wished to analyze three related dependent 
variables, requiring that the corrections obtained via 01 (or 
more correctly, SI) will not upset the relationship between the 
predicted values of the variables. Let the error in the predic- 
ted variables be denoted by Z(x,y), X(x,y), and Y(x,y), where 
(x,y) gives the spatial location and it is assumed that 

X(x,y) = k.Z (x,y) and Y(x,y) = k_Z (x,y). The subscripts x and 
X X ^ y 

y denote partial differentiation with respect to x and y 
respectively. Assume that the errors in the predicted values are 
stationary (that is, the statistics do not depend on (x,y)>, and 
have zero mean. Using EC.l to denote the expected value, or 
ensemble average, the spatial covariance function for Z, as a 
function of "lags" s and t, is 

R(s,t) = ECZ(x,y)Z(x+s,y+t) 3 * ECZ (x-s,y-t ) Z (x ,y) 3 . 

The latter equality follows from stationari ty. Under the assump- 
tion that the order of partial differentiation and the expected 
value can be interchanged, the cross covariance functions and the 
covariance functions for the derived variables are found in the 
manner illustrated here. ' Of course, it is assumed throughout 
this paper that the necessary derivatives exist. 
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ECZ <x ,y) X (x+s,y+t) 3 <■ ECZ (x ,y ) k <x+s,y+t) 3 = 

ECZ (X ,y> kjZg(x+s,y+t) 3 = k^ECZ <x ,y) Z (x+s,y+t) 3^ = kjR^(s,t) , 

whi 1 e 

ECX <x ,y) Y <x+s,y+t) 3 = ECX (x ,y) k^Z^ (x+s,y+t) 3 = 
ECX(x,y)k 2 Z^(x+s,y+t) 3 = k^ECX (x ,y) Z (x+s,y+t ) 3^ = 
k^ECkjZ^ (x-s,y-t) Z (x ,y) 3^ = k^EC-k^^Z^ (x-s,y-t) Z (x ,y) 3^. = 

-k. k_ECZ (x ,y) Z (x+s,y+t) 3 . = -k,k_R. (s,t) . 

Note that while the covariance functions are symmetric, the 
cross covariance functions are antisymmetric, which accounts for 
the sign change that comes from changing the order of the product 
in the expected value. This means, among other things, that the 
cross covariance must be zero at zero lag values. This behavior 
can be seen in the function plots in Bergman C113 and Scniatter, 
et. al. C103. 

2.2 Some Necessary Properties 

In order for the covariances of the derived functions and 

the cross covariance functions to exist, certain conditions must 

be satisfied by the function R(s,t). These have been alluded to 

by Buell C173, and are given in Julian and Thiebaux C183, where 

R (s) R (s) 

lim — is finite, and limC — - R (s) 3=0 

s>0 ® s»0 ® 

s in this equation represents the lag distance (lag in the above 
2 2 1/2 

was (s + t ) ), and R(s) is an isotropic covariance function. 

When one considers that R <0> must be zero, the first limit is 
the definition of the second derivative at s=0, hence existence 
of the limit means that the covariance function must be twice 
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differentiable at s=^0. The second limit then says that the 
second derivative is continuous at s^O. Thus the theorem given 
by Julian and Thi^baux can be simplified: 

Theorem 1: If R(s) is an isotropic covariance function for Z in 

two dimensions, then the covariance functions for the 
partial derivatives of Z exist at s=0 if and only if R(s) is 
twice continuously differentiable at s=0. 

2.3 Anisotropic Functions 

It has been contended that isotropic covariance functions do 
not adequately model the forecast error statistics and that gains 
can be made by using anisotropic functions. See Thi^baux 
C 133 , C 193 , C203 , and Thi^baux, et. al . C213 for development and 

discussion of product forms of covariance functions. Use of 
products of single dimensional functions has the advantage of 
carrying over desirable properties to higher dimensions, as well 
as being able to use essentially one dimensional structures and 
techniques. On the other hand, perusal of contour plots of 
product functions show that zero crossings of the functions occur 
along grid lines, and it is easy to see this will always happen. 
This may be undesirable behavior, and almost certainly it is not 
the kind of behavior seen in the error statistics. 

Another form of anisotropy is possible, one which results 
from scaling differently in two orthogonal directions, then using 
an isotropic functions in the scaled variables. This would 
result in the zero crossings in the contour plots of the function 
being ellipses with axes in the two directions, and all contours 
having the same shape. The eccentricity of the ellipse is a 
measure of the anisotropy of the error statistics. It would be 
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Basy to allow rotation along with the scaling to obtain ellipses 
of constant "distance" with any axis orientation. For a discus- 
sion of this type of anisotropic correlations, see Seaman C223, 
and Buell and Seaman C233. The properties of any such functions 
are those of isotropic functions, of course, since the anisotropy 
arises purely from a rotation and scaling. 

2.4 Isotropic Functions 

The use of isotropic functions in two or more dimensions 
that have been derived from one dimensional considerations can 
possibly lead to nonpositive definite functions. For example, 
Ripley 1241 (p 11) quotes a result of Mat6rn C25D , which gives a 
lower bound for isotropic positive definite functions in several 
dimensions. The result means that positive definite functions xn 
two dimensions are necessarily bounded below by -0.403 (the 

minimum value of Jq(s) ), while in three dimensions the bound is 
—0.218 . Thus any oscillatory positive definite function in one 
dimension that takes on values less than -0.403 cannot be an 
isotropic positive definite function in two dimensions. A posi- 
tive definite function with parameters to separately control the 
oscillation frequency and the decay can probably be made into a 
nonpositive definite isotropic function in two dimensions. For 
example, an exponentially damped cosine function, f (s) = cos (as) 
exp(— bs), can be made nonpositive definite by suitable choice of 
parameters, say a = 5 and b = .1 . This result also applies to 
other candidates for isotropic correlation function models, as 
will be shown later. 

There is a one-to-one correspondence between covariance 
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■functions in one dimension and isotropic covariance -functions in 
two dimensions. Using the so called "turning band” method « 
Matheron C261 gives a way o-f generating an isotropic d-dimen- 
sional covariance -function -from a one dimensional covariance 
■function. The relation is 



C . (s> 

U 



- 



, , ,, 2, <d-3> /2 . 

(vs) ( 1-v ) dv , 



where K is a constant that is unimportant -for our purposes. In 
two dimensions, this gives 



C2(.) 






, , ,, 2 , - 1/2 . 

(vs) ( 1-v ) dv . 



It is possible to invert Matheron 's relation to show a one-to-one 
relationship. A sketch o-f the inversion process follows. 



Employing a change of variables in the previous expression gives 

•s 



C„(s) = K / C. (t> (s^-t^> dt , 



1 



2 2 

then making further change of variables, s = x, and t * y» yields 






K / CCj(y‘ 

J r \ 



/ 2 , ,- 1 / 2 - , ,- 1/2 .. 

) (2y ) ] (y-x ) dt . 



1/2 - 1/2 

This is Abel's equation for K (y ) (2y) , and the solution 



is well known (see Hochstadt, C271), given by 



„ , 1/2, 1/2 d /_ , 1/2,, ,-1/2 . 

Cj (X > = K X ^ ^ dy , 



where K' is a different constant. Substitution for s and t once 



again, gives 



■i'=> - d7i=7/ 



2 2 - 1/2 

C2(t)(s -t ) '^2t dt . 



The correspondence between covariances in one dimension and 
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three dimensions is easier to invert, and is given by Ripley 
C243. There the relation is 

CjCs) =/Cj(vs) dv, and (s) = ^ CsC^(s)3 . 

•'o 

While this characterization o-f multidimensional isotropic 
covariance functions is interesting, and can in fact be used to 
generate isotropic multidimensional covariance functions, it does 
not easily answer the question as to whether or not a particular 
one dimensional function is an isotropic positive definite func- 
tion in more dimensions. One way to answer such a question is to 
use the characteri zation of positive definite functions as 
Fourier transforms of probability density functions (or alterna- 
tively, as functions whose Fourier transform is positive) . The 

Fourier transform of an isotropic function C(s) in two dimensions 

1 /2 

becomes (essentially) the Hankel transform of s C(s). It may 
be considerably easier to look at the one dimensional Fourier 
transform. It would then be useful to have a sufficient condi- 
tion on the Fourier transform of the function which would 
guarantee it is an isotropic positive definite function in two 
dimensions. Such a condition will now be derived. Let (s) be 
a positive definite function of one variable. From the charact- 
erization in the previous section, it can be shown 



(2.1) (s) = / cos(rs) h(r) dr. 



for some probability density function h(r) (i.e. , h(r)>0, with 
integral equal 1). The problem is then to determine the condi- 
tions that will make Cj(s) the two dimensional Fourier transform 
of an isotropic probability density function. Such a transform 



11 



is necessarily isotropic. A -function g<s) is sought so that 

.a> 

J^(rs) s g<s) ds 



< 2 . 2 ) 



Cj <r) 



/; 



This expression is inverted using the Hankel transform, giving 



(2.3) 



g(s) = f r (r) dr. 



Then, using (2.1) in (2.3), and interchanging the order of 
integration, followed by integration by parts yields 



g(s) = / J^(sr) r ( / cos(tr) h(t) dt)dr 

0 •'0 
• CD ^ <D 

h(t) ( / J^isr) r cos(tr) dr) dt = 



I h(t) sin(tr) dr) dt = 

0 •'o 

^ <D 

/ h ' (t) ( / Jq ( sr ) sin (tr ) dr ) dt , 



and then , 



(2.4) 



g (s) = 



/ 2 2 1/2 
- / h ' (t) / (t -s ) dt 



— 1/2 

The last equality uses the Hankel transform of r sin(tr). 

In order for g(s) to be a probability density function it 
must be nonnegative with integral equal to one. It is easy to 
show (again, interchanging the order of integration) that the 
integral is equal to one. It is more difficult to show necessary 
and sufficient conditions for the nonnegativity of g(s). The 
above relations summarized gives: 

Theorem 2: A sufficient condition for Cj (s) to be a valid iso- 

tropic covariance function in two dimensions is that h(t) be 
a monotone decreasing (h'(t)<0) function. 
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This condition sseems unnacasKAri 1 y rastrictiva, and dlf-fi- 
cult to usa since the condition is on tha Fourier cosine 
transform of <s) rather than (s) itself. Nonetheless, the 
condition can be used to show the following interesting results. 

I. Consider the exponentially damped cosine function, 

C(s) s= cos(as)e 

The Fourier cosine transform of this function is 

b (b^-<-a^+t = ) 



h (t) = F(C) (t) = 



Cb=+(a-t) =:Cb=+(a+t) =1 



2 2 

Inspection of h'<t) shows that if b > 3a , it is nonpositive for 
all t, and hence h(t) is monotone decreasing under that 
constrai nt . 

II. Consider the second order autoreqressi ve (SOAR) covar- 
iance function, 

C(s) = Ccos(as) + <b/a) si n <as) le 

The Fourier cosine transform of this function is 

2b (b=-4-a=) 



h(t) => F(C) (t) = 



Cb=+(t-a)=]Cb=+(t+a)=] . 



2 2 

Inspection of h'(t) reveals that if b > a , it is nonpositive 
for all t, and thus h (t) is monotone decreasing under that 
constraint. We see that each of the above C(s) is an isotropic 
positive definite function, hence is a covariance function if the 
appropriate inequality on the parameters is satisfied. 

III. Consider the special case of the damped cosine 
function 

C(s) = CA + (l-A)cos(as) :/ <1 + (bs)^)^^^ . 

The Fourier transform of this function is 
h(t) = F<C) (t) = 

(2b)~^f (1-A) CKq( It-al/b) + K^dt+al/b)! + AK^(t/b)> . 
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Because the modi-fied Bessel -function becomes unbounded as the 
argument tends to zero, -for A 1 , the Fourier trans-form must be 
increasing as t approaches a through values smaller than a. For 
t greater than a, and possibly -for some values smaller than a, 
the -function is decreasing. Thus the su-f-ficient conditions given 
above are not met, and it is easy to -find con-figurations o-f (x,y) 
points and parameter values A, a, and b for which the resulting 
"covariance" matrix is not positive definite. The two dimen- 
sional Fourier transform of C<s) <the Hankel transform of 
1/2 

s C(s) > has thus far gone unsolved, so it is presently unknown 
if there are parameter values (other than for A = 1) that will 
yield a positive definite function. 

IV. Consider the Bessel function The Fourier 

transform of this function is 



h(t) = F(C) (t) 



0 , t<a 

. 1 / 2,,. 2 2 , 1/2 .. 

t /<t-a) ,t>a. 



This function is easily seen to be monotone decreasing for t>a, 
and thus the Bessel function J^(as) is an isotropic covariance 
function in two dimensions. Application of this relation 
requires attention to some technical details because of the 
infinite Jump discontinuity at t=a. 

The above results concerning several functions proposed for 
use as isotropic covariance functions in two dimensions is 
useful. The lack of results and empirical evidence against the 
damped cosine being positive definite negate the results noted in 
the next section where we see that the fitting power of the 
function is very good. These aspects of the function will be 
discussed further in the next section. 
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2.5 Summary 

Th« contents of this section contain some useful information 
for the construction of isotropic positive definite function* and 
testing of functions for positive definiteness. When possible, 
the two dimensional Fourier transform of Cj (s> can be used to 
decide whether or not the function is positive definite. When 
the two dimensional Fourier transform cannot be obtained in 
closed form. Theorem 2 can give some information if the one 
dimensional Fourier transform is available in closed form. While 
the sufficient condition given by Theorem 2 is not necessary, it 
has been shown to be useful in investigating some functions which 
have been proposed for use as isotropic covariance functions in 
two dimensions. 

3.0 Some Experiments with Isotropic Covariance Functions 

3. 1 Background 

The work reported in this section is intended to help deter- 
mine something about the overall fitting properties of various 
suggested covariance functions. The term "overall fitting orop- 
erties" is meant to include not only the ability of the function 
to model a reasonably complicated true covariance function, but 
also its performance when used in a statistical interpolation 
scheme with several different observation patterns. 

The approach for this project was to begin with published 
data from an actual case, and then construct a covariance func- 
tion using a least squares fit to the data from a certain class 
of covariance functions. This model is used to define the 
"truth" model. Functions from other classes were then fit to the 
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same data, again in the least squares sense, and the performance 



of these "assumed" covariance functions were measured against 
that of the optimum model. The results to be discussed give some 
insight into what classes of functions have adequate fitting 
ability for modeling actual forecast error statistics, and also 
show how much skill is lost (in the idealized case) by use of 
inaccurate covariance functions. 

The results given here consist of representati ve plots of 
assumed correlation functions together with the correlation 
function defined as "truth", and contour plots of some of the 
resulting expected errors. The tables show expected root— mean- 
square (erms) errors (relative to the standard deviation of the 
background error) over three grids of points and associated 
observation locations. The expected errors were computed as in 
Seaman (1983). The results obtained with various assumed corre- 
lation functions in the SI scheme are discussed in detail. 
Additional plots are given in Franke C283. 

3.2 The Model Correlation Function 

The data for the covariance function was obtained (by hand) 
from Lonnberg C291. The data taken was plotted points from a 
covariance function of the type used by European Center for 
Medium-Range Weather Forecasts, in this case a five term (i.e., 
n»5) Bessel series of the form 



the radius of the region of interest. This function is positive 



n 






A.J^(s*k./R) + Art , 

f 1 0 1 0 ’ 




zero of the Bessel function J^(s), and R is 
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definite as an isotropic function in two dimensions provided the 
coefficients are ail positive. In Lonnberg , R was 2000 km. 

In this workf distance was measured in degrees, and the radius 
was scaled to 30°. 

Least squares fits to the data by functions of the type 
<3.1) for four, five, and six— terms were computed. While the 
original paper indicated a series with five terms generated the 
data, it was found that six— terms yielded all positive coeffi- 
cients and a significant reduction in the residual over five 
terms. Thus, it was decided to adopt the six term series as the 
"truth" covariance function. This six term series would also be 
marginally harder to approximate using other classes of 
covariance functions. The data and the fits using four and six 
terms are shown in Figure 1, and the coefficients are given in 
Table 1, along with other data. The intercept values of the 
approximations were 0.8270 and 0.8592 for four and six terms, 
respectively. This occurs because the data reoresents the 
spatial correlation of the background plus observation error, 
thus the interCBot is a function of the ratio of the standard 
deviations of background and observation error. The effects of 
this kind of discrepancy will be discussed in section 3.3. The 
correlation function for background error is the approximation 
normalized to have value one at s=0, of course. 

3.3 The Grid and Observation Point Sets 

Three grids and associated point sets were selected for 
studying the expected errors of statistical interpolation schemes 
based on various assumed covariance functions. All were based on 
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ths approximate locations o-f radiosonde data <-from Wahba and 
Wendelberger C303 and Ghil, et.al. C3i3) within the selected 
grid- Each grid covered a region that was 30° in longitude and 
20° in latitude, and the three were chosen to represent a dense 
observation set, a partially dense observation set, and a sparse 
observation set. The regions correspond to the middle United 
States with 36 observations, the eastern United States and west- 
ern Atlantic Ocean with 25 observations, and the middle Atlantic 
Ocean with 3 observations. For reference purposes, the three 
regions will be referred to as the MUS (Mid-US), EC (East Coast), 
and MA (Mid-Atlantic) regions. The regions and the observation 
locations can be seen in Figures 2—3, parts (b) , (c), and (d) , 

respectively. The regions were gridded at 2.5° intervals for 
purposes of computing expected errors, although the erms errors 
given in Table 1 are only over the interior grid points to mini- 
mize edge effects. Use of interior grid points for this purpose 
is valid since on a sphere it is not necessary to interpolate to 
the boundary points. For contouring purposes the fields were 
interpolated to finer grids using bicubic spline interpolation. 

3.4 The Assumed Correlation Functions and Results 

The families of assumed correlation functions fell into 5 
classess (i) Bessel function, (ii) negative squared 
exponential (sometimes called Gaussian), (iii) autoregressi ve of 
order two, (iv) autoregressi ve of order three, and (v) damped 
cosine. They will be discussed in turn, along with the results. 
Plots of the assumed correlation functions, along with the 
"truth" correlation function, are shown in part (a) of Figs. 2-3. 
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For fitting purposes, each included a multiplicative parameter 
that determined the s»0 intercept, and was subsequently dropped 
to obtain the correlation function. The value of this parameter 
is of interest, however, because dropping it shifts the curve 
(upward) to pass through the point (0,1), and thus different fits 
may be shifted by differing amounts, which ultimately affects the 
fit to the background error correlation function. 

Recall that the erms errors given in Table 1 are given as a 
fraction of standard deviation of the background error. The 
ratio of the standard deviation of the observation errors to the 
standard deviation of the background errors was 1/3. 

<i) Bessel Function 

The reference expected errors were computed using the actual 
correlation function model, given by Eq. (3.1) with coefficients 
as given in Table 1. The results are given in Table 1, and are 
the smallest expected errors that can be obtained using a correc- 
tion to background scheme, that is, they are truly optimum. The 
correlation function is shown in part (a) of Fig. 2, while the 
contour plots of the expected error for each of the three grid/ 
observation sets is shown in parts <b), (c), and (d) . 

The results using a four term Bessel function are given in 
Table 1. Because the intercept of the fit to the data is 0.8270 
versus 0.8592, normalization to value one produces a curve that 
is predominately above that for the model correlation function, 
especially for small distances. The result of the poor approxi- 
mation for small distances is most pronounced over MUS. The 
effect was small over the sparse part of EC and over MA. 
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(ii) Negative Squared Exponential (NSE) 

The NSE has been recognized as inadequate for modeling error 
covariances -for some time, and the results obtained here confirm 
that. The assumed form of the function was 

(3.2) A + (l-A)e"^®''^^'* . 

This function is positive definite as an isotropic function in 
two dimensions for 0<A<1. 

The initial fit was not obtained by least squares, but 
simply by attempting to fit the model correlation reasonably well 
for small distances, taking A=0. The fit is reasonably good up 
to about 6°, and quite poor at greater distances. The expected 
errors are similar in magnitude to the expected errors for the 
four term Bessel function, except over MA, where the errors are 
larger. However, since the errors over MA tend to be large 
anyway, the relative effect is not as great as one might expect. 

The second attempt was by least squares for the parameters A 

and b. Because the NSE is too flat near the origin, this process 

yielded an intercept value of 0.8060, shifting the correlation 
function so that it is entirely above the model correlation 
curve. This results in even poorer performance over MUS and EC 
than the previous model, due to the inaccurate representati on for 
small distances. The performance over'MA was better than the 
above. 

Due to the poor performance (compared to the above) obtained 
by adding a constant to the basic NSE it was decided to attempt 

to find a better fit by trial and error. No claim is made about 

any optimality for this function. The results in Table 1 demon- 
strate that it is probably not possible to obtain good results 
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ovsrall with a function from tho NSE family, and certainly not 
for the present model correlation function. 

(iii) Autoregressi ve. Order Two (SOAR) 

The SOAR model has been suggested as appropriate by Yudin 
C321 , Thi^baux C133 and this is supported by simulations due to 
Balgovind, at.al. C333. This is the model that is being incor- 
porated into the U. S. Navy NWP models. The formula given here 
includes a constant term which is not part of the SOAR model, but 
which has been noted to improve performance considerably 
(Thidbaux, et.al., C213), and those results are confirmed here. 
The SOAR function with additive constant is 

(3.3) A + ( 1— A) Ceos (as) + (b/a) sin (as) 3e 
This function is positive definite (in two dimensions) whenever 
a< b , and 0<A<1. In ail cases investigated here, and as has been 
reported elsewhere, (e.g., Thiebaux, er.al., C213), the parameter 
a tends to be essentially zero. In this case the function re- 
duces to 

(3.3a) A + (l-A)Cl + bsle”^^® . 

The initial attempt was a least squares fit to che data with 
A=0. The intercept obtained was 0.7977, with the resulting 
correlation curve then being considerably above the model corre- 
lation curve between 0° and 15°. The performance was only 
slightly better than with any of the previous correlation 
functions. It was then decided to attempt a least squares fit 
with the intercept constrained to be 0.8592, the same as. obtained 
for the model correlation function, but again with A=0. Table 1 
shows marginal improvement for all three grid/observati on 
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patterns. A third attempt included A in the least squares -fit, 
with no constraint. This resulted in a much closer match to the 
model correlation -function, although the Intercept o-f 0.8441 
moved the assumed correlation curve above the model curve -for 
much o-f the interval. The -fit and resulting expected error 
contours are shown in Figure 3. Table 1 shows there is consid- 
erable improvement over all previous results, the most improve- 
ment being -for MUS, and the least -for MA. 

<iv) Autoregressi ve, Third Order <T0AR) 

The use o-f the TOAR model has been investigated by Thiebaux, 
st al . C213, including an additive constant. The formula is 

(3.4) A + ( 1-A) C (acos (as) + bsin(as))e + ce , 

where the coefficients a, b, and c are functions of a, b, and c 
given by 

a = (3b'^-a'^-c‘^) ac/D , b = (b‘^-3a"^-c‘^) bc/D , 
c = -2(b^+a^)ab/D , where D = (3b^-a^-c^) ac-2 (b^+a^) ab . 

It is unknown what restrictions (beyond 0<A<1) on the parameters 
are required to ensure the function is positive definite as an 
isotropic function in two dimensions. 

The data was fit by least squares with the TOAR function 
(3.4). The intercept was 0.8651, which resulted in the curve 
being slightly below the model correlation curve over most of the 
range. Overall, the fit was quite close and better than any of 
the previously discussed functions. The results in Table 1 show 
very close agreement with the optimum possible for all three of 
the grid /observation sets. 
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(v) Damped Cosine 

The damped cosine -function has been suggested by Thidbaux 
C191 and Seaman and Hutchinson C343. The formula is 
<3. 5) CA + (l-A)cos(as) 5/Cl + (bs)^]*^ . 

It is unknown whether the function is positive definite as an 
isotropic function in two dimensions, but the evidence in section 
2.4 (for c=0.5), while inconclusive, seems to indicate it is not. 
In practice, of course, the function may be positive definite 
when the observation points are restricted to certain regions. 

The data was fit with the function (3.5), under the restriction 
c=0.5. The intercept was 0.SS65, which resulted in a very slight 
raising of the curve relative to the model correlation function. 
The resulting fit is excellent for small distances and very good 
over the entire range. Table 1 shows that this function gives 
the best results of all the functions tested. 

(vi) Variations 

The expected error computations for a number of variations 
of the above functions were also performed. The principal 
variation was to fit the data only over the first half of the 
interval, (0°,15°). The effect of this was to generally (though 
not always) increase the erms errors over MUS and EC, while not 
affecting the results over MA. In the damped cosine, the 
exponent c was chosen by least squares, along with the other 
parameters, and resulted in a slightly better fit to the correla- 
tion function, especially at larger distances. However, the 
coefficient A was slightly greater than 1. Whatever the positive 
definiteness properties of the function, having A>1 will certain- 
ly make it non— positive definite. Although no graphical results 
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are shown, the coefficients and erms errors are given in Table 1 
for the additional assumed correlation functions. 

3.S Sensitivity of the SOAR Model to Parameter Mi sspecifi cation 

In order to determine more completely the characteristics of 
the SOAR model, some additional calculations were made to deter — 
mine the effect of mi sspecifi cation of the parameters in the 
correlation function or the ratio of the standard deviations of 
the observed and background error. The results can be summed up 
rather quickly: The scheme is mostly insensitive to such 

variations. Figure 4 shows a family of 4 correlation functions, 
#4 being the SOAR plus constant discussed in the section 3.4, 
with the others having smaller correlations at a given distance. 
Figure 5 shows the expected RMS error for each of the four as the 
"assumed" correlation, when the "true" correlation function is 
#4. With the exception of the sparse MA grid, the expected 
errors are relatively stable under significant perturbations. 
Figure 6 shows the sensitivity to the assumed error ratio, and 
once again, it is observed that the expected errors are quite 
stable. 

4.0 Conclusions 

The principal conclusion to be drawn is that the correlation 
family used in practical analysis should embody a sufficient 
number of parameters to fit the forecast error statistics reason- 
ably well. Further, it is most important that the data be fit 
accurately for small distances. In order to ensure a better fit 
for small distances, it may be worthwhile to enforce the inter- 
cept of the correlation for the background plus observation 
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errors ii_ the ratio of standard deviations of the two errors is 



known accuratel y . The effect of scaling to obtain the correla- 
tion function, and the apparent shift up or down can possibly be 
compensated for by artificially varying the ratio of background 
to observation errors, as well, although it seems more desirable 
to enforce this ratio in the correlation function fitting 
process. 

As noted above, clearly the most important region for the 
fit to the correlation function to be accurate is for small 
distances. Over the sparsely observed region, MA, and to a 
lesser extent over the EC region, the overall erms errors were 
only slightly affected by the assumed correlation function. In 
the case of the MA region it is noted that the error contours are 
relatively unaffected except near the observations. Since the 
errors in the remote part of the region dominate the overall 
error, the choice of assumed correlation function has relatively 
small influence. On the other hand, over the densely observed 
region, an accurate fit at small distances was most important. 

The NSE correlation function, while not performing well, illus- 
trates the above nicely. For the first NSE entry in Table 1, 
even though the fit is poor for distances of more than 6°, the 
erms errors over MUS and EC are smaller than the best fit (last 
NSE entry) due to the more accurate fit for small distances by 
the former function. Of course the arms errors over MA are 
poorer for the first case due to the very bad fit at large 
distances. 

There appear to be several good candidates for use as two 
dimensional isotropic correlation functions, including SOAR, 
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TOAR, and damped cosine, given by Eqs. (3.3), (3.4), and (3.5), 

respectively. While the -fitting power for the latter two are 
greater (there are a greater number o-f parameters -for those two) , 
the choice o-f SOAR seems reasonable and adequate -for a number o-f 
reasonss (1) The SOAR (with the additive constant) embodies a 
su-f-ficient number o-f parameters to allow oscillation and decay 
with distance. (2) The SOAR has some credibility as the spatial 
correlation -function o-f an innovation process. However, results 
are -for one dimension rather than two, except -for the results 
cited previously in Balgovind, et.al. C331. (3) The SOAR was 

demonstrated here to be positive de-finite as an isotropic 
-function in two dimensions, under a mild restriction on the 
parameters. <4) While the TOAR is also the spatial correlation 
-function (again in one dimension) -for an innovation process, 
based on this limited study it does not appear to be signi-fi- 
cantly better than SOAR. (5) The positive de-f ini teness 
properties o-f the TOAR are not known, although it is certainly 
positive de-finite as an isotropic function in two dimensions 
under some restrictions on the parameters. (6) Although the 
fitting ability of the damped cosine seems to be at least as good 
as the TOAR, and it is positive definite in one dimension, evi- 
dence indicates it may not be positive definite as an isotropic 
function in two dimensions, regardless of parameter restrictions. 
The availability of other acceptable alternatives seems to make 
it prudent to preclude the use of the damped cosine in practical 
si tuati ons. 

It is pointed out that all of the functions except the four 
term Bessel function and the NSE perform very well. Table 1 
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shows, for example, that the SOAR is only a little more than 17. 
of the standard deviation of the background error poorer than 
optimal over hUS and EC, and less than .IX poorer over HA. 

Finally, it is noted that within the SOAR family, SI is 
quite insensitive to misspecif ication of the correlation parame- 
ters, even to an extent such that the correspondence would appear 
to be much less between two members of the family than between it 
and a fit by the NSE. Thus it could as important to choose the 
correct family of correlation functions as well as to model 
properly within that family. In addition, misspecif ication of 
the ratio of standard deviations of the background and observa- 
tion errors has a rather small effect on the skill of the method. 

This work has focused only on the univariate problem, 
whereas in practice such schemes are applied to the multivariate 
one. Further work is necessary to determine whether the nice 
results obtained here carry over to the multivariate case. A 
further investigation of the effect of wind observations on the 
analysis of pressure height and wind fields is anticipated. 
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Parameter Values and arms Errors 
Model Correlation is Six Term Bessel 



Assumed 

Correlation Parameters erms 



Function 


a,b,c 


A 9 


MUS 


EC 


MA 


6 T Bessel 




0.2474 
0.3335 
0. 1844 
0. 1031 
0.0362 
0.0554 
0.0400 


0.2667 


0.3752 


0.7483 


4 T Bessel 




0.2811 

0.3090 

0.2213 

0.0930 

0.0956 


0.3046 


0.4088 


0.7503 


NSE 


10.0 


0.0 


0.3047 


0.4184 


0.7822 


NSE 


14.88 


0.3200 


0.3688 


0.5282 


0.7631 


NSE 


10.0 


0.2500 


0.3098 


0.4158 


0.7541 


SOAR 


0.0 

0. 1215 


0.0 


0.3034 


0.4022 


0.7562 


SOAR 


0.0 

0. 1374 


o 

a 

o 


0.2931 


0.3968 


0.7562 


SOAR 


0.0 

0.2055 


0.2722 


0.2780 


0.3859 


0.7491 


TOAR 


0.4732 

0.3828 

0.0914 


0. 1974 


0.2717 


0.3794 


0.7485 


Dmpd Cos 


0.4749 
0. 1367 
0.5000 


0.9592 


0.2686 


0.3779 


0.7486 


NSE„ 


15.0 


0.0 


0.3619 


0.4414 


0.7649 


NSE 


12.31 


0.3205 


0.3474 


0.4299 


0. 7593 


SOAR 


0.0 

0,2654 


0.3758 


0.2743 


0.3825 


0. 7495 


TOAR* 


0.4468 
0. 1482 
0.0052 


-5.9965 


0.2697 


0.3801 


0. 7514 


Dmpd Cos* 


1 . 2236 
0. 1507 
0. 5000 


1 . 0027 


0.2749 


0.3734 


0.7491 


Dmpd Cos 


0. 7009 
0.2069 
0.3753 


1.0105 


0.2692 


0.3779 


0.7484 


Dmpd Cos* 


0.7987 

0.2350 

0.3317 


1.0147 


0.2706 


0.3784 


0.7485 



Table 1 

These correlation -functions were obtained by least squares 
■fit over the interval <0°,15°) 
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FIGURE CAPTIONS 



Figure 1: The data points and least squares fits by -four and six 

term Bessel -functions. 

Figure 2: (a) Six term Bessel correlation function (true and 

assumed). (b) Expected root-mean-square error contours for the MUS 
grid and observation point set for the correlation functions shown 
in (a). (c) As in (b) for the EC grid and observation point set. 

(d) As in (b) for the NA grid and observation point set. 

Figure 3: (a) Six term Bessel correlation function (true) and 

second order autoregressi ve plus constant correlation function 
(assumed) . (b) , (c) , and (d) As in corresponding parts of Fig. 2 

for the correlation functions shown in (a). 

Figure 4: Four 2nd order autoregressi ve correlation functions, 
as in Eq. 3.3a, with (b,A) values: #1 - (0.5, 0.0); #2 - 
(0.3, 0.0); #3 - (0.3»0.15); #4 - (0.2055,0.2722). 

Figure 5: Expected RMS errors when the "true’’ correlation is 

function #4 with various assumed correlation functions for each 
of the three grid/observation sets. 

Figure 6: Expected RMS errors when the assumed ratio of the 

standard deviations of the observed to background error is 
varied. Actual error ratio is 1/3. 
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FIGURE 1 
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